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Abstract 



Neuronal growth cones are the most sensitive amongst eukaryotic cells in responding to directional 
chemical cues. Although a dynamic microtubule cytoskeleton has been shown to be essential for 
growth cone turning, the precise nature of coupling of the spatial cue with microtubule polariza- 
tion is less understood. Here we present a computational model of microtubule polarization in 
a turning neuronal growth cone (GC). We explore the limits of directional cues in modifying the 
spatial polarization of microtubules by testing the role of microtubule dynamics, gradients of reg- 
ulators and retrograde forces along filopodia. We analyze the steady state and transition behavior 
of microtubules on being presented with a directional stimulus. The model makes novel predic- 
tions about the minimal angular spread of the chemical signal at the growth cone and the fastest 
polarization times. A regulatory reaction-diffusion network based on the cyclic phosphorylation- 
dephosphorylation of a regulator predicts that the receptor signal magnitude can generate the 
maximal polarization of microtubules and not feedback loops or amplifications in the network. Us- 
ing both the phenomenological and network models we have demonstrated some of the physical 
limits within which the MT polarization system works in turning neuron. 

Key words: Microtubule-dynamics; neuron; growth-cone; polarization; reaction-diffusion; gra- 
dient; 



Introduction 



A long-standing goal in developmental biology has been to understand the process of neuronal 
path finding that results in the complex wiring seen in the nervous system. Studies have unravelled 
multiple cues that provide directional signals to the axon. Once the directional cue is detected by 
a the growth cone (GC) at the end of an axon, the intracellular cytoskeleton undergoes polarized 
growth. The signaling networks that transduce an extracellular directional cue into cytoskeletal 
polarization during a neuronal GC guidance show striking similarities with those in highly migratory 
cell types such as slime molds Dictyostelium, neutrophils and migratory fibroblasts (j5"5"|) . This 
diversity of cell types suggests an evolutionarily conserved set of mechanisms in cell guidance. 
Experimental work in nerve GCs has shown both actin and microtubule dynamics to be essential 
(j45p . While many studies have investigated the role of actin in filopodial and lamellipodial dynamics 
in neuronal growth cones, the precise role of MTs is less well studied (reviewed in (j23j) ) . MTs 
have been hypothesized to function like a cellular compass, which by random fluctuations searches 
directional space. In the presence of a cue, a reinforcement of the direction by transport of actin 
polymerization factors and membrane vesicles leads to turning of the GC ([1]). Thus MT polarization 
appears to be an important early event in GC turning, but consists of complex interplay of chemical 
and mechanical effects. 

The chemical regulation of MTs in neuronal growth cones is through modulation of MT dynamic 
instability- the property of MTs to transition from growing state to a shrinking state (catastrophe) 
and vice-a- versa (rescue) (|28l S3J). This dynamics has been shown experimentally to be essential 
for growth cone turning (|T3l l4"5l |62HM|) and is most sensitive to the frequency of catastrophes ([S7|) . 
Important regulators of length dynamics is a family of neuron specific proteins of the stathmin 
(opl8) family- stathmin 2 (SCG10) and stathmin 3 (SCLIP) (56). The protein SCG10 destabilizes 
microtubules by increasing catastrophe events and on phosphorylation results in stabilized MTs, 
corresponding to decreased catastrophe events (|24[ 158ft . in a manner similar to stathmin. SCG10 
levels in primary cultured neurons appear to be critical for maintaining a balanced level of MT 
dynamics to enable neurite extension (j47p and vivo experiments indicate the N-methyl D-aspartic 
acid (NMDA) receptor (NMDAR) activation regulates SCG10 phosphorlyation through ERK2 (|ig|) . 

Mechanical force generated by the retrograde flow of actin flowing inwards from filopodial tips 
regulates the number and velocity of MT ingress into filopodia in neuronal growth cones (|60|) . The 
coupling of the actomyosin driven flows with the MT system is thought to occur through discrete 
sites of binding by factors with either both actin- and microtubule-binding activity (e.g. formins) or 
complexes of factors that have such activity (12ip . An additional force acts in the opposite direction 
of the retrograde flow driven by cytoplasmic dynein which works to both translocate microtubules 
towards the neuronal periphery, and to hold them in place in the filopodia (|50l) . 

Theoretical models and simulations have been used to model some sub-cellular aspects of neu- 
rons providing insights into their dynamics. The cytoskeletal dynamics in growth cone polarization 
was modeled in filopodial dynamics were modeled to be randomly searching space (|65p . A model 
of microtubule assembly and tubulin transport could reconcile the constant elongation rates of 
neurites with tubulin assembly processes (66). Complex two-dimensional models which include the 
retrograde flow of actin suggest a non-random mode of MT invasion of filopodia (|27|) . A recent 
detailed model of retrograde actin flow and myosin contractility has revealed the balance of forces 
driving the flow in neurons (|15|) . However all these models ignore the chemical regulatory ele- 
ments that play a crucial role in neuronal cytoskeletal dynamics. Independently, chemical spatial 
gradients have been modeled in MT spindle assembly (|33|) and cell polarization (|54|) . Previous 
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modeling of reaction-diffusion gradients and their regulation of MT dynamics has been used to ex- 
plain the non-random growth of MTs in mitosis (4, 69). More recently, models of signaling proteins 
in neuritogenesis have demonstrated the role of feedbacks (J2j) and of cell shape in axon determina- 
tion (53). However a neuronal growth cone model that includes both regulatory reaction-diffusion 
components and the explicit MT dynamics and mechanics has not yet been developed. Given the 
accumulating evidence for the importance of microtubule regulation in growth cone guidance, a 
spatial model that includes both MT and actin interactions, combined with spatial dynamics of 
signaling proteins that modulate MT growth will be useful both from an experimental as well as 
theoretical perspective. 

Here we develop a model based on Aplysia growth cone geometry and dynamics with exper- 
imentally measured parameters. Such a model provides on the one hand a platform for testing 
experimental hypotheses that can lead to extremely involved experiments. We develop two kinds 
of models: a phenomenological MT stabilization model as well as an explicit biochemically realis- 
tic reaction-diffusion model of MTs in turning GCs. Some of the typical scenarios of GC turning 
by a diffusible cue and perturbations of the actomyosin system show good agreement with exist- 
ing experiments. Additionally we predict the minimal spatial and temporal detection limit. The 
maximal MT polarization angle shows good agreement with the maximal GC turning angle seen 
in experiments. The detailed biochemically realistic regulatory reaction-diffusion network model 
biologically is based on the cyclic phosphorylation-dephosphorylation of a microtubule destabilizer 
of the stathmin family, found in neurons. It predicts the most effective point of regulation in such a 
network. We have used both the phenomenological and network models to understand the physical 
limits within which the MT polarization system works in neuron that have detected a polarizing 
cue. 



Theory 

Growth cone and axon shaft geometry 

The neuronal growth cone and terminal axon shaft were modelled as a 2D objects. The geometry 
was inferred image analysis of published images of Aplysia bag cell neuronal growth cones (|12|) 
using Image J ()61j) and a routine developed in- house using MATLAB R2008b (Mathworks Inc., 
USA). A circle of radius 25/im is used to represent the growth cone based on dimensions of Aplysia 
growth cones. The growth cone is divided into concentric circles with the outer peripheral P- 
domain (radius 10 fJ,m), the intermediate region transition T-zone (radius 5 \xm) and the innermost 
central C-domain (radius 10 //to). The confining boundary is modeled as set of connected points 
with stiffness set to 100 pN/^m, making it effectively rigid. A rectangular region representing the 
axon-shaft emerges from the lower end of the circular growth cone. Based on reports of microtubule 



bundling by actomyosin forces in the growth cone neck (|11|) continuous spatial force field (Fig. 1A ) 
was implemented. For every fiber model-point lying in this field, the value of the force at the 
corresponding position in the field is added to the net force on the point. 

Microtubule orientation and dynamics 

Microtubules were implemented as elastic, non-extensible, polar fibers represented by a set of dis- 
crete model-points. The flexural rigidity (or bending stiffness) of microtubule fibers was set at 20 
pN-/irn 2 (22 , 142p . which corresponds to a persistence length of 4.67 mm. Bending of these fibers at 
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a shorter length scale is however possible based on the rigidity and force applied (29, 38j). To main- 
tain the simplicity of the model, we do not model microtubule breakage. The fibers were modelled 
to be dynamic based on the four-parameter model (|67|) : growth (v g ), and shrinkage velocities (v s ), 
frequency of rescue (fres) and catastrophe (feat)- The parameters were spatially inhomogeneous 
with dynamic instability parameter values different for each of the domains as described in Table 
[TJ These parameters affect only that plus-tip that enters it. An additional cortical catastrophe was 
modelled within lfj,m of the membrane periphery of the growth cone. 



A fixed number of free microtubules (Nmts = 30) are initialized with lengths randomly assigned 
between 16-24 [im by a uniform random number generator with a total angle of 106° symmetrically 
about the vertical axis in an arc centered below the center of the circular growth cone (Fig. 1A). 
This geometry is based on Aplysia growth cone micrographs. The modeled MTs are dynamic at 
both ends. The parameters of dynamic instability: frequencies of catastrophe and rescue, and 
velocities of growth and shrinkage are distributed inhomogeneously as a function of radial position 
in the growth cone (Fig. 1A). Inside the axon-shaft, minus ends of the modelled MTs undergo 
rapid catastrophe. The field of dynamic instability is set to a high f ca t value (de-stabilizing) in the 
P-domain, while the T-zone and C-domain have a stabilizing effect (Fig. 1A). 

The minus ends of the MTs were initialized 20 fim below the point of contact of the axon shaft 
with the growth cone in the axon shaft. The minus-ends are dynamic and undergo catastrophes if 
they move beyond 20/im radial distance from the center due to a spatial change in he catastrophe 
frequency to f ca t = Is" 1 , while the other parameters remain constant (shrinkage velocity v s = 
0.5fim/sec, rescue frequency f res = Is -1 and v g = 0/um/s). 



Force due to actin retrograde flow in the P-domain 

The actomyosin retrograde flow is modeled by discrete immobilized motors arranged radially along 



filopodial bundles that push MTs inwards towards the C-domain (Fig. IB) based on reports in 
literature of coupling-molecules that can bind to MTs and actin simultaneously. In filopodia the 
actin itself is flowing inwards in a myosin-dependent retrograde flow, hence generating many point 
forces that pull MTs into the center of the growth cone. Recent reports on the role of kinesins 
in MT growth and polarization in the P-domain suggest that both the actomyosin system, and 
a kinesin based retrograde force are involved in moving MTs inwards from the periphery (|51 j) . 
Assuming both these systems act along the same direction, the force vectors due to actomyosin and 
kinesin activity can be added up and resulting in a high stall force value (f stall) (Table [2]). Control 
calculations showed this force is sufficient to move MTs of a fixed length inwards with a retrograde 
flow speed that agrees with experimental measurements (60). 

Motors-MT interactions are modeled as a Hookean springs where force exerted by a motor 
on a point of the fiber is given by f ex = k * 5r, where k is the stiffness and 5r is the distance 
separating the motor and the point on the fiber. Motors were characterized by their positions, 
attachment/detachment rates (r attach A 'detach), speed of movement (v mo t), direction of movement, 
stall force (f stall) and stiffness (k) of the induced links. Detachment rate depends on the force gen- 

erated as modelled previously (|36|) as r detach = r detach * ex P fstaLl based on Kramers theory (i37|) . 
The filopodial arrays (n= 61) are arranged between and 180° with respect to the positive X-axis 
through the semi-circle of the GC with 30,000 motor complexes randomly distributed (uniformly) 
amongst the filopodial bundles ~ 500 motors per bundle). The speed of these motor complexes 
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was set to 0.07 fim/s while attachment-detachment rates were optimized (see Table [2]) to provide 
net translocations velocities in agreement with measurements in Aplysia growth cones (I60p . 

In the T-zone an analog of immobilized non- motile dynein motors is modelled with v a ttach — 5, 
^detach = 0.1 and v-mot = 0. These motors hold MTs but do not move the microtubules when bound 



(Fig. IB). 



Reaction-diffusion network model of receptor driven intracellular polarization 

A simple reaction-diffusion network consisting of three species: two cytosolic diffusive species as- 
sumed to be stathmin S and its phosphorylated form S* and a membrane immobilized receptor 
species (R). An optional positive- feedback interaction by growing microtubule plus tips (N t i ps ) and 
a weight of the feedback (u>i) was implemented. This feedback increases the S— > S* conversion 
in presence of growing plus-tips. S is hypothesized to be phosphorylated to S* , similar to stathmin 
phosphorylation, near the membrane by the receptor R, mimicking the NMDA-Erk2 pathway. The 
backward reaction converts S* to S homogeneously in the cytosol governed by a first order rate 
constant k 2 . In the positive- feedback model, growing MT plus tips can also enhance this forward 
reaction. The forward reaction is modelled using Michaelis-Menten kinetics as follows: 

d{S]_ h[R][Sr N k fb - Wl -[Sr . m v2[<?1 m 

~dT " ~K M1 + [5]"> - NUPS • K M2 + [S\<» + k2[S ] + Ds ' V [S] (1) 

The backward reaction is modelled using first-order kinetics as: 

d[S*} _ Mg][gr i N k ff ■ m • [S]" 12 u rc*l, D V 2r^ m 

~ k mi + [ S ]^ + Nups • k M2 + [5F " h[s ] + Ds * • v [s] (2) 



The terms n\ and n 2 correspond to cooperative Hill-coefficients for the forward reaction and 
the MT-dependent feedback terms respectively. These terms represent cooperativity in phospho- 
rylation (ni) or the positive feedback (r^), due to either the extracellular regulated kinase (ERK) 
pathway (39) or multiple sites of phosphorylation in stathmin (j40|) . Unless explicitly stated, the 
Hill- coefficient parameters were set to n\ = n 2 = 1. Diffusion coefficients of the two diffusing 
cytosolic species S and S* are given by D$ and D$* while the concentration gradient is indicated 
by the 2D Laplacian V 2 . 

The role of S in modulating MT dynamics was implemented by increasing the local catastrophe 
frequency f cat dependent on [S] using a phenomenological expression described previously based on 
in vitro data ([3]): 

feat = feat ■ exp(w 2 • ([5] - [Sua])) (3) 

Here f ca t is the default value (Table [I]), when [S] = [S]tot, [S]tot is the sum concentration of modified 
and unmodified stathmin, and w 2 = 0.3144, a scaling factor obtained by fitting to experimental 
data (3). The parameters used for the reaction-diffusion network are taken from literature and 
described in Table [SlJ 
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Methods 



Simulation of the microtubule model 

Our model is implemented and simulations conducted using Cytosim- a cytoskeleton simulation 
engine implemented in C++ (I52p . All simulations were run for 1200 seconds (or more). Typical 
simulations include 30 dynamic MTs, 30000 surface-immobilized motors in P-domain and 1000 
grafted motors in T-zone. The reaction diffusion system was simulated with spatial resolution of 
1/im. Positions of all objects and occurrence of events like movement, binding/unbinding, change 
in length of MTs (due to growth/shrinkage), change in MT end state (due to catastrophe/rescue) 
is calculated for all objects at each time step from the corresponding rate parameters. Simulations 
were run on a 12-core Dell T5500 workstation with 2.40 GHz processors using a publicly available 
parallel processing code (Distributed Parallel Processing Shell Script vers 2.85). Each typical run 
takes ~ 1500 seconds to complete. 

For each iteration total MT numbers, mean direction of MT distribution and circular spread 
of MTs were calculated in the P-domain and angular trends were analyzed in 12 angular sectors, 



15° each (Fig. 2C). The mean angle of MT distribution in a growth cone at a time point (t) was 
calculated as: 

«.<*) = 5^ (4) 

^tips 

where 6i is the mid-angle and N£ ips are the number of plus-tips in the ith sector. The corre- 
sponding circular standard deviation (a c (t)) is given by 



<Jc(t) 



UNI 

tips 



Time averages of mean angle {{9 av )) and circular standard deviation ((<r c )) were used as mea- 
sures to compare with experiments. In experiments where a directional bias was introduced, the 
difference between mean angles averaged at steady state before and after the bias (A(6 av )) were 
evaluated as a measure of MT polarization. 



Simulation of the reaction-diffusion model 

The 2D system of PDEs was solved on a finite lattice grid using by an explicit forward difference 
scheme (|35|) . The spatial interval Ax = \\xm and time interval At = O.Olsec were optimized by 



comparing simulation results with analytical solutions of known initial condition cases (Fig. S4) 



(j!6p . Zero- flux boundary conditions were used at the growth cone and axon-shaft boundaries 
(Fig. 1A). All preliminary calculations were performed in Octave (I20p and then implemented in 
C++. 



Analysis of experimental MT distributions 

Images of fluorescently labelled microtubules in Aplysia growth cones were analyzed in house using 
a MATLAB script (Mathworks Inc., USA) and ImageJ (NIH, USA) combined with manual counting 
of plus tips. 
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Results 



Microtubule angular distribution is uniform and fluctuating 



A two-dimensional model of the Aplysia growth cone with dynamic MTs has been developed to 
understand the quantitative details of spatial remodeling of MTs in the earliest events relating to 
neuronal growth cone turning. The model integrates MT regulatory and mechanical components 
as described in the model section. The simulation geometry with spatial dynamic instability and 



filopodial geometry is shown in Fig. 1 A A radial inward-directed retrograde flow and dynein based 
stationary coupling of MTs is modeled using a motor model geometry schematically illustrated 
in Fig. IB A simulation run of such a model shows MTs entering the P-domain in ~ 10 min to 



attain a uniform but fluctuating distribution (Fig. 2A, Movie SI). We also observe some MT 
buckling mostly in the T-zone. This geometry and dynamics are qualitatively similar to in vitro 
data of Aplysia growth cones. The quantitative dynamics of MTs entering the P-domain also show 
rapid steady state where on an average 2 /3 rd of the MTs remain in the P-domain (Fig. 2B). The 



fluctuations suggest that MTs can constantly search the P-domain for directional cues, primed for 
turning. The number of MTs in P-domain for two different values of initialized MT numbers (30 



and 60) show between 25-50 MTs on an average (Fig. 2C). This is in agreement with previous 
experimental observations I60|) . The spatial patterns of these MTs were analyzed by sampling 
the number of plus-tips in discrete angular periods for each time point (Fig. |2D[ ). The weighted 
mean angle for the MTs in the P-domain of a GC was calculated based on Equation [4] and plotted 
against time (Fig. 2E). After initial fluctuations reaches 6 av reaches a steady state value of ~ 90° 
in about 300 s of simulations, while the corresponding circular standard deviation (<r c ~ 7r/6) 
calculated from Equation [5] achieves a steady state value of ~ 45°. The time- aver aged mean MT 
angle ({0 av }) is as expected ~ 90° ((<t c ) ~ 30°), independent of the number of MTs simulated 



(Fig. 2F). Both the mean angle and the circular standard deviation are in good agreement with 
experimental measures (|114 l60|) . Simulations of MT density dynamics in four equal sectors of the 



growth cone for initial MTs N=30 (Fig. |3Aj), N=60 (Fig. (3Bj) and N=90 (Fig. [3C]) are all of the 
same order of magnitude as experiments with Aplysia growth cones (|32p . 

In order to test if the standard model of motor detachment based on Kramers theory affects our 
results, a simple model with constant rates of motor detachment in the P-domain was tested. The 
values rdetach = Is or rdetach = 5s _1 are the two extreme values from a two-state load dependent 
detachment model which was shown to better fit experimental data for kinesin (|19p . The mean MT 



number in the P-domain and the mean angle remain unchanged (Fig. SI), suggesting little effect 
of changing detachment rates. 

Thus our model simulations show qualitatively and quantitatively expected behavior in agree- 
ment with the available spatial and dynamic data. In the following sections we perturb this model 
system and compare it with experiments. 



A decreased retrograde flow rate results in increased MT flux in the P-domain 

The role of the actomyosin retrograde flow was explored in our modeled growth cones by analyzing 
its effect on the flux of MTs entering the P-domain. A reduction in the retrograde flow rate 
(from 4.2fj,m/min to 2.1/im/mm) by 50% results in a two-fold increase in MT-tip flux only in the 
peripheral 25% of the modeled growth cone (Fig. |3D[ ). Such a dramatic effect is not observed if 
larger proportions (50% and 75%) of the P-domain were sampled. This suggests the fluctuations in 
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the most distal region of the growth cone are most sensitive to such treatment. On the other hand a 
two-fold increase in retrograde flow rate results in an ~ 85% decrease in flux of MT tips in the distal 
25% of the P-domain. This effect is also seen in experimental analyses of MT dynamics which report 
MT fluxes only in the distal 25% of the growth cone. The effect of changing retrograde flow rates 
on absolute number of MTs in the P-domain shows small but statistically significant differences 



(Fig. 3E). The relative proportion of simulated MTs (N=30) compared to those entering the P- 
domain also progressively diminishes if we observe the most peripheral 25% of the growth cone. 
Experimental measurements of the flux of MT plus-tips in with neuronal growth cones perturbed 
by blebbistatin - a myosin II ATPase inhibitor- indeed showed a 50% decrease in retrograde flow 
rates results in a two-fold increase in the flux of MTs entering the P-domain in the peripheral 25% 
of the P-domain 0120 . This correspondence with experimental data shows our model of retrograde 
flow can be used for further studies on MT dynamics modulation. 

A spatial catastrophe frequency change is most effective in changing the polar- 
ization angle of MTs 

Recent evidence has been accumulating for the role of MT polymerization modulation in growth 
cone turning. We have systematically tested the role of MT dynamics in MT polarization, an early 
even in growth cone turning. Values of dynamic instability were modified in only one half of the 
P-domain of the simulated growth cone. The result of this spatial bias was evaluated by plotting a 
time course of the 9 av {t). Although fluctuations in the mean MT-angle still remain large, a small 



change in mean angle is observed from the unbiased base (Fig. 4A). Both the actual simulation 



(Movie S2 ) and a time averaged angular distribution of the MTs before and after bias in the control 



case compared to the biased case represent this change more graphically (Fig. 4B). In order to test 
the parameter to which the system was most sensitive, the bias was systematically introduced in 
each of the parameters of dynamic instability (feat an d fres)- The change in mean angle (A{8 av )) 
of MTs in the outer half of the P-domain was evaluated and compared to the control (cat/res) 



(Fig. 4C ) . Increasing catastrophe frequency turns the mean angle of MT distribution away from 
the direction of bias, while the increase in rescue frequency turns the MT distribution towards 
the bias. The ten- fold increase in f res , ten-fold decrease in f ca t and ten- fold increase in f ca t show 
the significant changes in mean angles compared to the unbiased case (Two-way Student's t-test 
p < 0.05). A ten-fold decrease in f res does not result in a significant change in mean angle of 
MTs. To further examine the minimal range within which the rescue and catastrophe frequencies 
change mean angles of MT, step-wise increases in f res and f ca t were made. Statistically significant 
chang es in A(B av } were observed even with two- and four-fold changes in fcat and fres values 



(Fig. 4D). These results suggest even a small change in dynamic instability can bias a the MTs in 
the given geometry, however changes in the fcat values produce the most dramatic effect on A(9 av ). 
These results are consistent with previous non-spatial models of MT dynamics that demonstrated 
theoretically as well as experimentally (in mitosis) that the mean MT length is most sensitive to 
changes in f ca t (|67|) . This also sets the stage for testing the spatial limits of MT polarization in 
our model growth cones. 

Spatial and temporal limits to MT polarization sensitivity 

In order to test the sensitivity of the MT polarization system in the modelled growth we modeled 
different extents of the polarizing cue and evaluated the output. The sensitivity was estimated in 
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two ways- (i) the change in mean angle (A{8 av )) as a function of the angular extent of the bias, 
and (ii) the time taken for the adaptation, i.e. the new steady state. In simulation f ca t values 
were decreased or f res values were increased ten-fold in sectors of the P-domain. These sectors 
were all centered around 120° and the angular extent was varied between — 60°. The change in 
mean angle A{9 av ) was reduce from ~ 15° to ~ 5° in response to both changes in catastrophe and 
rescue frequencies. The change in mean angle with decreasing spread of f res showed a seemingly 
steady polarization angle of > 10° for bias extents even as small as 30°. A rescue frequency bias 
less than 30° also resulted in a small but statistically significant changes (Student's t-test) in the 



mean angle (Fig. 5A). On the other hand reducing the spatial extent of f ca t resulted in a steady 



reduction in A(9 av ) with ever decreasing spatial extent (Fig. 5B). This result appears to contradict 
the observed maximal sensitivity of MT length to the catastrophe frequency. It is likely to be a 
result of the increased catastrophe frequencies at the GC periphery, which are better compensated 
by a higher rescue frequency. This prediction can also be experimentally tested, and would suggest 
alternative sensitivity of the MT system in realistic cell geometries as compared to in vitro or 
mitotic measurements. Additionally it suggests that a measurable MT polarization cannot occur 
if the bias is smaller than 20° in angular extent for catastrophe- and 10° for rescue-frequencies. 

In order to test the temporal sensitivity of the polarization of MTs, one half of a growth cone 
was biased by changing either ca t or res abruptly. The time taken for adaptation was estimated by 
fitting the dynamics of the mean angle of MTs in the growth cone {0 av {t)) to a simple saturation 



kinetic expression (Fig. 5C): 



,(t) = e av (0) + A0 av (t) ■ (6) 



where av {t) is the mean direction at time-step t, 8 av (0) the mean direction before the bias is 
introduced, A{0 av ) is the net change in direction after the bias, n is the Hill coefficient, and Kt 
is the half-time within which polarization is half complete. The change in Kt from the time of 
imposition of the bias is AKt. A smaller value of AKt indicates a faster response to the stimulus. 
Ten-fold increases and decreases (independently) in catastrophe and rescue frequencies produced 
AKt values ranging from 30-45 s. A ten-fold decrease in rescue frequency could not be fitted by the 
routine. This appears to indicate a characteristic time for half-growth cone bias experiments for MT 
polarization to be in the range of ~ 30s. This could be experimentally tested in a carefully designed 
experiment. At the same time, the phenomenological gradient of stabilization of MT growth in 
space must have an underlying mechanistic basis. Hence we developed a spatially explicit reaction- 
diffusion model of a potential MT regulatory protein in the framework of the previously developed 
spatial MT dynamics in the neuronal growth cone. 

A cyclical phosphorylation-dephosphorylation reaction-diffusion network leads 
to MT polarization 

In order to test the role of a biologically realistic factors that can change microtubule dynamics, a 
reaction-diffusion spatial model of an MT regulator was developed, based on the role of the protein 
SCG10, the stathmin (op 18) homolog. We have simplified the network of extracellular regulatory 
kinase (ERK) dependent phosphorylation of stathmin to consider its core element, namely MT 
spatial regulation. In our model stathmin (S) phosphorlyation (S*) occurs at the membrane due 



9 



to a receptor (R) while dephosphorylation of can occur throughout the growth cone, resulting 
in a gradient of f ca t (Fig. 6A). An additional feedback term (u>i) regulates further amplifies the 
phosphorylation reaction. Both S and S* are modeled to be diffusible through the growth cone 
geometry, while R is a parameter in the model, as seen in Equations [I] and [2] In order to test 
numerical stability of our integration scheme, the diffusion of a single species was compared to 
that in an analytical model (Fig. S4). The modulation of MT dynamics occurs in our model by 
modifying the f ca t values based on Equation [3j The resultant steady state distribution in such a 
simple model of phosphorylation-dephosphorlylation cycles demonstrates the differences in spatial 
extent of S and S* in response to a stimulus in one half of the growth cone as seen along the radial 
(Fig. 6B ) and central (Fig. 6C ) axes. The time course of gradient formation by depletion of S in 



one half of the GC is visible in Movie 

In order to test which component of this reaction scheme might affect the polarization of 
MTs maximally, we compared the A(8 av ) for receptor signal strength and nature of the reac- 
tion kinetics. We tested the effect of changing the phosphorylation reaction of S to S* transition 
from a simple first-order (KMi/[Stot] = 2), intermediate (KMi/[Stot] = 1) or zero-order saturated 
(KMi/[Stot\ = 0.1) scenario. The effect of changing receptor concentrations [R] (1, 10, 100, 200 
/iM) was simultaneously tested for each of the KMi/[Stot] values. We find a high receptor signal 
values lead to most significant changes in A(6 av ) (Fig. |7A[ ). This is also borne out by the mean 
steady state post-bias angular distributions of MTs (Fig. |7Bj ). All reactions were carried out with 
the Hill-Coefficient in the phosphorylation reaction set to n% = 1 (Equation [I]) . For the situation 
of a high receptor stimulus ([R] = 100/iM), when n\ was set to higher values (2, 4 and 10) the 
polarization angle A(9 av ) did not change appreciably (Fig. [S3] ) , indicating a negligible effect of the 
cooperative term. 

The feedback from MTs to promote the phosphorylation reaction was simulated and compared 
to the model without feedback, in terms of the final readout for the f ca t distributions. The simple 
model demonstrates a slight but clear directional bias, while the feedback model in contrast shows 
local inhomogeneities in f ca t, corresponding to the presence of growing plus-tips (indicated by 
arrows) (Fig. 7C). As a result the mean post-bias MT distribution is visibly polarized in the simple 
reaction case as compared to the strong feedback model (with both positive-feedback w\ = 200 and 
Hill-coefficient m = 10) (Fig.[7D|). 

Taken together, these results suggest a model of a stathmin-like phosphorlyation-dephosphorylation 
mechanism without feedbacks can indeed polarize microtubules in a neuronal growth cone geometry. 
It remains to be seen if such a gradient and downstream MT effect can be measured experimentally. 



Discussion 

Cytoskeletal dynamics in GC turning has been increasingly seen as critical for understanding neu- 
ronal GC guidance ()17p since multiple signaling pathways converge on the actin-microtubule system 
(|18p . However, very few of the details of this coupling between the signaling mechanisms with the 
cytoskeleton in neuronal growth cones are known. Previous studies have focussed either only on 
the mechanics of cytoskeleton or the dynamics of the signaling system. Our study for the first 
time explicitly considers this coupling between microtubule mechanics and signaling gradients in 
a neuronal growth cone. We focus on the role of the spatial regulation of MTs in a turning GC. 
Another feature of our model is its comparison with previously published data of MT polarization 
in an Aplysia neuronal growth cone, thus making many aspects of our model biologically relevant 
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and testable. All input parameters of microtubule dynamics and actin retrograde flows are taken 
from neuronal studies, except in the case of microtubule-motor mechanics and stathmin, where the 
parameters are derived from in vitro studies. 

Our neuronal MT-polarization model is based upon the following assumptions: (1) a static 
extracellular guidance cue, (2) single step phosphorylation of the MT regulatory protein, ignoring 
intermediate steps (3) a phenomenological non-linear mapping of stathmin concentration and mi- 
crotubule dynamics and (4) the two-dimensional geometry of the growth cone with the mechanics of 
retrograde pushing of MTs. During the process of neuronal GC turning the assumption of a static 
guidance cue holds true since the time scale of remodeling of the MTs is an order of magnitude 
slower than receptor polarization. For instance MT polarization in Aplysia growth cones occurs 
within ~ 10 min (121} I32D . while receptor polarization of NMDAR takes only seconds (25, [68]) . We 
use this difference in time-scales to focus on the events slower than membrane-receptor polariza- 
tion, i.e. downstream signaling and cytoskeletal regulation. The downstream signaling from these 
receptor signals either through GTPases or through kinase-phosphatase systems have a conserved 
modular structure, permitting a simplification of the cascade as a simple two-component cyclical 
reaction (j^ |34"|). Using such a relatively simple single-step reaction cycle, we can also test multi- 
ple modes of amplification and feedback in a manner that is tractable. Such amplification steps 
are thought to be part of the mitogen activated kinase (MAPK) extracellular-regulated kinase 
(ERK) pathwav (|39p . upstream of the cyclical reaction. Although the Hill-cooperativity modeled 
in Equations [T] and [2] has not been measured in neuronal GCs, this model allows us to explore the 
multiple levels at which potential amplification of the reaction network can occur. The modeled 
two-dimensional geometry was assumed to be a reasonable approximation of the growth cone based 
on micrographs of Aplysia bag neuron growth cones which are typically 20-30 jjam in diameter but 
only 2-3 \im in height in the P-domain (32). 

Using this model we compare our simulations to previously published experimental data. It is 
clearly seen that the mean angle and spread of the modeled neuronal MTs are strikingly similar 
to the experimentally measured values from Aplysia growth cones in the absence of a directional 
cue [60). The simulated quadrant- wise time dynamic of MT densities in the peripheral 25% 
were plotted for GCs with the initial number of MTs as 30 (Fig. [3A]), 60 (Fig. [tlB]) or 90 (Fig. 



3C). Agreement with experimental measures of MT density in the peripheral P-domain (32) and 



quadrant-wise dynamics (62) is seen only for the MT numbers 60 and 90. However the effective 



MT-flux rates in simulation (Fig. 3D , Table S2 ) are similar compared to experimental values (|12j) . 
The comparison of the change seen in MT flux in the peripheral 25% of the P-domain after a 50% 
reduction in retrograde flow rates is further evidence that the model of microtubule dynamics and 
geometry appears to be valid, given the available data. The converse experiment of the increase in 
flux rates has been performed by siRNA inhibition of dynein in the GC which oppose the retrograde 
flow. In experiments this results in an ~ 8 fold reduction in MT's in the peripheral P-domain 
(|50p , a difference not seen in our model. However the interpretation of the experiment as well as 
differences between rat superior cervical ganglia neurons from Aplysia neuron sizes, could partly 
cause the difference. We also find a spatial bias in the retrograde flows velocities can influence MT 
polarization in either two- fold increase or decrease in retrograde velocities (Fig. S2), making our 
results comparable to selective inhibition of kinesin in only one half of the GC which resulted in 
growth cone turning (|51j) . 

A surprising property to emerge from our model simulations was the buckling of microtubules 
in the T-zone (Fig. 2 A ) . It results from the combination of the growth cone geometry, immobilized 
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motor complexes of the T-zone, and actin retrograde flow in the P-domain. Such MT buckling 
has been previously observed in experimental measurements (60), and is thought to also result in 
breakage of MT fibers, which then either polymerize or are translocated towards the P-domain. In 
previous iterations of the simulation, it was found that the T-zone localization of immobile motor 
complexes was critical to maintain the orientation and spread of MTs as they enter the P-domain. 
Hence our model would appear to suggest a functional role for such a set of complexes, that at first 
sight appear passive. Changing these complexes into minus-end directed active motors increases 
the buckling of the MTs. A further quantification of this buckling both in experiment and theory 
might prove useful. These molecular motor assemblies might also play a role in the proposed 
polarity sorting of short microtubules in neurons proposed in qualiitative models before by Baas et 
al. (J5j[6|). Our model could be used in future to test the specifics of such processes. 

The sensitivity of the MT polarization system in response to changes in dynamic instability 
parameters was measured to give an idea of the ideal points of regulation. We found the system to 
be most sensitive to the frequency of catastrophe. However statistically significant changes in mean 
MT angle were also observed as a result of fold-changes in both catastrophe and rescue frequencies. 
The change in mean MT angle is observed to be 20 — 30° for rescue and catastrophe frequencies. 
The experimental mean turning angle observed for Xenopus embryonic (|10p and spinal neurons 
(j43p as well as rat spinal neurons (jSJ) were also in the range of ~ 20°. We consider this to be a 
strong, but qualitative validation of our model of MT polarization. 

Regarding the spatial sensitivity of the system, the model predicts that the attractive or repul- 
sive cue needs to be spread across minimally 30° angle of the growth cone. Most experiments of 
GC turning involve a pipette assay where the chemical that provides the directional cue is applied 
extracellularly centered 45° to the GC and angular spread appears quite large (~ 90°) and is dif- 
ficult to control in such vitro assays (l57|) . However bead based assays (j4"6|) suggest an alternative 
method for estimating the role of spatial extent of the cue. The larger question of the in vivo 
scenario might also further complicate matters with multiple cues simultaneously stimulating a 
single GC. We cannot resolve this in our simulation. Our model does however predict a lower limit 
to to the MT polarization system in the neuronal GCs, which can be experimentally tested. The 
guidance of neurons was shown to increase linearly with signal to noise ratio (SNR) both in theory 
and experiment with rat dorsal root ganglion (DRG) explants (@9j). Our results are qualitatively 
comparable with previous data of increased MT polarization with increasing receptor concentration 
and broadening of the receptor signal. However a more detailed receptor-ligand dynamics model 
will be required for a direct comparison between the two models. 

We further resolve our spatial MT polarization model of GCs by including a stathmin-like reg- 
ulator of MT dynamics and its cyclic phosphorylation-dephosphorylation, since stathmin has been 
shown to modulate MTs and affect neurite outgrowth in rat hippocampal neurons (|47p . The map- 
ping from phospho-stathmin concentration to MT dynamics parameters however has been measured 
only in Xenopus oocytes ([3|) and we use the empirical fit to this data as a start point. This biologi- 
cally realistic scenario is used to test the effect of different network topologies such as a feedback from 
MT tips to the reaction network and saturated reaction kinetics. We find these amplifications do not 
improve the angle of turning, as compared to the simple cyclic phosphorylation-dephosphorylation 
reaction system. The feedback from MTs intact appears in our system to amplify noise and cause 
'spurious' polarization even in the absence of a cue. On the other hand we do find that receptor 
activation strength proves to be a much better means of strengthening the polarization of the MT 
system. This result reveals an important design principle of such a system- namely amplification 



12 



steps appear more effect at the receptor signaling level. This is consistent with recent findings in 
a simple model of MT feedback based GABA receptor polarization that amplifies receptor based 
signaling (7). Our study goes beyond this to analyze the explicit regulation of MT dynamics by 
both chemical and physical mechanisms. 

In conclusion, the model presented here predicts various mechanical and chemical properties of 
MT cytoskeleton regulation in neuronal GC turning. Validation with existing data demonstrates 
the utility of this model and points the way for additional experimental testing. The predicted 
physical limits to MT polarization in GCs is important, since it reveals potential design principles 
of neuronal GC turning. Additionally recent developments in microfluidics and micro patterning 
techniques <\31\ |59|) suggest means of experimental testing of these predictions. The recent in vivo 
evidence from neuronal regeneration by MT stabilization (I26|) shows the growing need to understand 
the role of microtubule polarization at a system level in neurons. We believe this study is a step in 
this direction. 
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Tables 



Location of MT plus-end 


P-domain 


T- and C-domain 


Reference 


v g (ptm — s^ 1 ) 


11.5 • 10~ 2 


6.4 • 10~ 2 




v s (fim — s -1 ) 


20 • 10~ 2 


5.1 • 10~ 2 




feat _1 ) 


1.7- 10~ 2 


0.5 • 10~ 2 


•}•> 


fres (s ) 


2.3 • 10" 2 


9.5 • 10" 2 


5) 



Parameters of Dynamic Instability 



Table 1: The dynamic instability parameter values were chosen for the P- and C- 
domains and T-zone based on published values from Aplysia neuronal GCs. 



Parameter 


P-domain 


T-zone 


Reference 


No. of motors 


30,000 


1,000 


optimized 


Speed v mo t (fim/min) 


4.2 







Binding rate r attach (s _1 ) 


5 


5 


optimized 


Unbinding rate r detach (s^ 1 ) 


1 


1 


optimized 


Maximum stalling force f stall 


15 


15 


4X3HX1I5X) 


(pN) (Kinesin and actin- 








retrograde flow) 








Attach distance (/jm) 


0.1 


0.1 


optimized 


Motor stiffness k (pN/fxm 


100 


100 


(US), d 



Motor parameters 



Table 2: The parameters of surface-immobilized motors were chosen for the P-domain 
and T-zone to fit the retrograde flow rates and effective translocation speeds of MTs 
reported (60). 
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Figure Legends 

Figure [TJ 

Geometry of growth cone MT-motor system (A) A semi-circular geometry represents the 
growth cone, which is divided into three concentric regions: P-domain, T-zone and C-domain. The 
C-domain is connected to an axon shaft. A MT stabilization zone in the P-domain promotes MT 
polymerization. Immobilized plus-end directed motors are arranged in sixty-one symmetric radial 
arrays spaced at intervals of 3° in the P-domain. Non-motile motors similarly immobilized are 
localized in the T-zone. The reaction-diffusion system is bounded in the light and dark regions 
(a subset of simulation space). An effective force field in the GC-neck and axon shaft represents 
a compressive force (white arrows). The spatially distributed dynamic instability parameters v g , 
v s , feat and f res for the plus- and minus-ends of are plotted along the GC-axon axis. (B) The 
geometry of the modeled MT-motor interactions in the T-zone with non-motile complexes and in 
the P-domain with plus-end directed motors is schematically represented. 

Figure [2j 

MT dynamics in an unbiased growth cone (A) The time-series of MTs entering the P-domain 
(the region with higher catastrophe frequencies) in simulation snapshots at 0, 600 and 12000 s is 
represented in 2D plots and (B) as a time course of MTs entering the P-domain and achieving a 
fluctuating steady-state. (C) At steady-state the number of MTs in the P-domain in simulations 
with 30 or 60 initialized MTs was compared to images from experiments (n=2) I60j) . (D) The 
geometry of the growth cone angular sectors is schematically represented. (E) This is used to 
calculate the time-dependent mean angle of MTs (9 av (t))and its standard deviation (a c (t)). (F) 
The simulated mean angle and spread at steady-state (t > 400s) was compared to experimental 
distributions (n=2) (f6"U|) . 

Figure [3} 

Effects of perturbing the number of MTs and retrograde translocation rates MT den- 
sities in the peripheral 25% of the growth cone are plotted against time for initial nucleated MTs 
numbering (A) 30, (B) 60 and (C) 90. (D) The number of MTs entering the P-domain per second 
show a 50% change when standard retrograde translocation speeds (4.2 (im/min) are changed by 
50%. This effect is visible in the most-distal quarter of the P-domain and not as pronounced when 
larger elements (50 — 75%) of the P-domain are sampled. (E) The absolute number of MTs in the 
P-domain are less dramatically affected by changes in retrograde translocation speeds. Error bars 
are standard deviation (Iterations n=10). Stars indicate statistically significant differences from 
the control case. 

Figure [4} 

MT polarization in a biased growth cone (A) The time-dependent mean angle (0 av (t)) of 
MTs shows a change in mean direction when a directional bias is introduced in the left half of 
the GC at 800s by decreasing f cat ten-fold. (B) The angular frequency distribution of the mean 
number of MTs at steady state shows a clear polarization in the biased case. The sectors are 15° 
apart while the radial distance from the origin indicates the number of MTs (Iterations n=10). (C) 
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The mean change in angle of MTs after the bias (A(9 av )) is the greatest for a a ten- fold increase 
in the f ca t values in one half of the GC, but the change is also statistically significant (star) for 
an increase in f res and decrease in f ca f. (D) Fold changes in f ca t and f res biases (Iterations n=10) 
show a statistically significant changes (star) for two-fold or greater changes in both parameters. 
Error bars indicate standard errors (n=10). 

Figure [5j 

Spatial and temporal sensitivity Decreasing angular sectors of the GC were biased by changing 
values ten- fold of (A) f res (increased) and (B) f cat (decreased). The sectors were centered around 
120° (striped) with respect to the GC semi-circle. Less broad biases in f cat resulted in a steep 
reduction in MT polarization, while it remained fairly similar for f res biases as narrow as 10°. 
Error bars indicate standard errors (n=10). (C) The time taken for the half-maximal change in 
mean MT angle after a bias is introduced (AKt) was estimated by fitting the time-course of mean 
MT angle (9 av (t))- (D) The most rapid change in MT mean angle was produced by a ten-fold 
reduction in f ca t, while a similar reduction in f res produced almost no change. 

Figure [6j 

Spatial reaction network (A) A schematic view of the phosphorylation of the regulatory protein 
S by a membrane receptor signal R converts it to S* . S itself increases the frequency of catastrophe. 
The dotted lines indicate the additional feedback of MTs that promotes the S to S* transition, 
further stabilizing MTs. The intensity of brightness (dark: low, bright: high) in the GC plots 
the spatial extent of S. (B) Lines indicating the (C) radial and (D) vertical (GC to axon) axes 
along which concentrations were plotted. Solving the simple cyclic reaction network, the initial and 
steady-state concentration profiles of S and S* show a gradient from the GC periphery inwards 
along both (C) radial and (D) vertical axes. The corresponding f ca t gradient profiles result from 
the concentration gradient of the network. 

Figure [7} 

The amplification of the signaling system (A) The polarization of MTs in response to zero- 
order (KMi/Stot = 0.1), intermediate (KMi/Stot = 1) and first-order (KMi/Stot = 0.1) reaction 
kinetics are compared to increasing values of the receptor signal (n=30). (B) The schematic polar 
frequency distribution of the number of MTs shows the greatest polarization in the MTs in response 
to a first-order reaction kinetic with the highest receptor signal. (C) Spatial plots of S distribution in 
the GC are compared for two networks- the simple cyclic reaction without feedback (R = 200//M, 
KMi/Stot = 2, n\ = 1) and the growing MT-tip driven positive feedback loop (with additional 
parameters w\ = 200, KM2/Stot = 0.1, 112 = 10). (D) The polar plots at steady state of number 
of MTs display a highly polarized distribution for the simple cyclic reaction without feedback as 
compared to the strong feedback network. 
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Supplementary Tables 



Pararr 


eEteascription 


Value 


Unit 


Reference 


[S)tot 


Total stathmin con- 
centration 


8 


liM 


© 


D s 


Diffusion coefficient 
of S 


15 






D s * 


Diffusion coefficient 
of S* 


15 


[im 2 s~ l 


HMD 


[R] 


Receptor concentra- 
tion 


[1, 10, 100, 200] 


\iM 




W\ 


Weight of positive 
feedback 


[1, 10, 100] 


\iM 




W 2 


Weight of f cat depen- 
dence on [S] 


0.3144 






h 


Forward rate con- 
stant 


0.1 








Backward rate con- 
stant 


0.1 


s- 1 


- 




Feedback rate con- 
stant 


0.1 






Kmi 


Michaelis-Menten 
constant (Km) for 
the forward reaction 


[0.8, 8, 16] 


fj,M 


_ 


Km2 


Michaelis-Menten 
constant (Km) 
for the feedback 
reaction 


8 






m 


Hil Coefficient of the 
forward reaction 


[1, 2, 4, 10] 






ri2 


Hill Coefficient of the 
feedback reaction 


[1, 10] 







Reaction diffusion network parameters 



Table SI: The parameters of the reaction diffusion scheme were chosen either based on 
previous experimental measurements or order of magnitude estimates. 
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Region of P- 


Circumference 


MTs crossing 


MTs entering the 


domain con- 


(fim) 


into the P-domain 


P-domain per 10 


sidered (w.r.t 




per minute 


fim — min~ l 


periphery) 








100% (r=15 fim) 


47.1 


14.4 


3.06 


75% (r=17.5 fim) 


54.95 


12 


2.18 


50% (r=20 fim) 


62.8 


6.6 


1.05 


25% (r=22.5 fim) 


70.65 


3.6 


0.51 



Flux of MTs entering the P-domain 



Table S2: The parameters of the reaction diffusion scheme were chosen either based on 
previous experimental measurements or order of magnitude estimates. 
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